6.6.1bioclim
## 代码地址：
https://jcoliver.github.io/learn-r/011-species-distribution-models.html#additional-resources
##
R中物种分布模型的简要介绍
杰夫·奥利弗
2020年7月20日
借助R套件套件，从纬度和经度坐标预测物种范围变得越来越容易。本入门教程将向您展示如何将坐标数据转换为范围图。
学习目标
安装用于物种分布建模的软件包
使用bioclim方法运行物种分布模型
在地图上可视化模型预测
物种分布模型正在成为越来越重要的工具，以了解生物体如何应对当前和未来的环境变化。物种分布模型（SDM）的方法和文献数量在不断增长，因此鼓励您查看“其他资源”部分中的一些资源。dismo包装的装饰图案特别有用，Jeremy Yoder的介绍是另一个很好的起点。在本教程中，我们将使用公开可用的数据来构建，评估和可视化仙人掌仙人掌的分布模型。
入门
在执行任何操作之前，我们将需要设置工作区，下载示例数据并安装运行模型和可视化其输出所必需的其他软件包。
工作区组织
因此，首先，在工作区中创建一对文件夹：
dir.create(path = "data")
dir.create(path = "output")
优良作法是保持输入（即数据）和输出分开。此外，任何存放在output文件夹中的工作都应该完全是一次性的。也就是说，数据和我们编写的代码的组合应允许我们（或其他任何人）再现任何输出。
示例数据
我们正在使用的数据是对仙人掌（Carnegiea gigantea）仙人掌的观察。我们使用的是全球生物多样性信息机构GBIF的部分记录。您可以从https://tinyurl.com/saguaro-obs下载数据；将其保存在data您在上述步骤中创建的文件夹中。
安装其他R软件包
接下来，将需要安装另外五个R软件包：
install.packages("dismo")
install.packages("maptools")
install.packages("rgdal")
install.packages("raster")
install.packages("sp")
模型的组成部分
物种分布模型背后的基本思想是采用两种信息源来模拟预期发生某种物种的条件。两种信息来源是：
发生数据：通常是在其中已观察到感兴趣物种的纬度和经度地理坐标。这些被称为“状态”数据。一些模型还利用“缺失”数据，这是已知没有发生该物种的位置的地理坐标。缺少数据比较难获得，但是某些建模方法需要这些数据。在本课程中，我们将使用您先前下载的saguaro的发生数据。
环境数据：这些是环境的描述，可以包括温度和降水的非生物测量以及生物因素，例如是否存在其他物种（例如掠食者，竞争者或食物来源）。在本课程中，我们将重点介绍WorldClim提供的19种非生物变量。与其从WorldClim下载数据，不如从dismo包中使用函数来下载这些数据（见下文）。
数据和质量控制
我们将通过加载所需的这五个库来启动脚本。当然，还要在脚本的最上方添加一些信息，说明脚本的作用以及负责人！
# Species distribution modeling for saguaro
# Jeff Oliver
# jcoliver@email.arizona.edu
# 2018-02-27
library("sp")
library("raster")
library("maptools")
library("rgdal")
library("dismo")
您很有可能已经看到一些红色消息显示在屏幕上，特别是在加载maptools或rgdal库时。这是正常现象，只要所有消息中都不包含“ ERROR ”，您就可以在这些消息中嗡嗡作响。如果加载库确实导致错误消息，请检查以确保库已正确安装。
现在我们已经加载了那些软件包，我们可以使用以下getData功能下载生物气候变量数据：
bioclim.data <- getData(name = "worldclim",
var = "bio",
res = 2.5,
path = "data/")
我们提供了getData四个关键信息：
name = "worldclim"：这表示我们要下载的数据集的名称
var = "bio"：这getData表明我们要下载所有19个生物气候变量，而不是单独的温度或降水量测量值
res = 2.5：这是我们要下载的数据的分辨率；在这种情况下，度数为2.5分钟。对于其他解决方案，可以通过?getData在控制台中键入来检查文档。
path = "data/"：最后，这将设置文件下载的位置。在我们的例子中，这是data我们在开始时创建的文件夹。
另请注意，将文件下载到data文件夹后，它们将被读入内存并存储在名为bioclim.data
# Read in saguaro observations
obs.data <- read.csv(file = "data/Carnegiea-gigantea-GBIF.csv")
# Check the data to make sure it loaded correctly
summary(obs.data)
## gbifid latitude longitude
## Min. :2.021e+08 Min. :26.78 Min. :-114.0
## 1st Qu.:1.453e+09 1st Qu.:32.17 1st Qu.:-111.4
## Median :1.571e+09 Median :32.28 Median :-111.1
## Mean :1.567e+09 Mean :32.16 Mean :-111.3
## 3rd Qu.:1.677e+09 3rd Qu.:32.38 3rd Qu.:-111.0
## Max. :1.806e+09 Max. :34.80 Max. :-109.3
## NA's :3 NA's :3
请注意，和列中有三个NA值。这些记录对我们没有任何用处，因此我们可以将其从数据框中删除：latitudelongitude
# Notice NAs - drop them before proceeding
obs.data <- obs.data[!is.na(obs.data$latitude), ]
# Make sure those NA's went away
summary(obs.data)
## gbifid latitude longitude
## Min. :8.910e+08 Min. :26.78 Min. :-114.0
## 1st Qu.:1.453e+09 1st Qu.:32.17 1st Qu.:-111.4
## Median :1.571e+09 Median :32.28 Median :-111.1
## Mean :1.575e+09 Mean :32.16 Mean :-111.3
## 3rd Qu.:1.677e+09 3rd Qu.:32.38 3rd Qu.:-111.0
## Max. :1.806e+09 Max. :34.80 Max. :-109.3
现在，当我们查看obs.data数据框时，没有NA值，因此可以继续进行操作。
为了使物种分布模型更加简化，了解物种在地理上的分布范围是很有用的。我们将找到一般的纬度和纵向边界，并存储此信息以备后用：
# Determine geographic extent of our data
max.lat <- ceiling(max(obs.data$latitude))
min.lat <- floor(min(obs.data$latitude))
max.lon <- ceiling(max(obs.data$longitude))
min.lon <- floor(min(obs.data$longitude))
geographic.extent <- extent(x = c(min.lon, max.lon, min.lat, max.lat))
在进行任何建模之前，最好通过在地图上绘制点来对发生数据进行真实性检查。
# Load the data to use for our base map
data(wrld_simpl)
# Plot the base map
plot(wrld_simpl,
xlim = c(min.lon, max.lon),
ylim = c(min.lat, max.lat),
axes = TRUE,
col = "grey95")
# Add the points for individual observation
points(x = obs.data$longitude,
y = obs.data$latitude,
col = "olivedrab",
pch = 20,
cex = 0.75)
# And draw a little box around the graph
box()
看起来不错！
建立模型并可视化结果
现在我们的出现数据看起来不错，我们可以使用生物气候变量来创建模型。我们要做的第一件事是将考虑范围限制在合理的地理区域内。也就是说，出于我们的目的，我们并不是要在全球范围内模拟仙人掌栖息地的适宜性，而应在西南地区进行建模。因此，我们可以将生物数据变量限制在发生数据的地理范围内：
# Crop bioclim data to geographic extent of saguaro
bioclim.data <- crop(x = bioclim.data, y = geographic.extent)
# Build species distribution model
bc.model <- bioclim(x = bioclim.data, p = obs.data)
## Error in .xyValues(x, as.matrix(y), ...): xy should have 2 columns only.
## Found these dimensions: 400, 3
哦哦 这不好。看来我们传递给的数据bioclim格式不正确。提示出现在错误消息的第二行：中## Found these dimensions: 400, 3。这是指obs.data数据帧，它确实具有400行和三列。从以下文档中bioclim（请参阅?bioclim控制台中的内容）：
用法
bioclim（x，p，…）
参数
x Raster *对象或矩阵
p两列矩阵或SpatialPoints *对象
…其他参数
因此，无论我们传递给什么，p都应该只有两列。让我们对其进行修改，obs.data使其只有两列。第一列是GBIF标识符，我们不需要它，因此我们使用负运算符（即减号）将其删除。然后我们可以运行物种分布模型。
# Drop unused column
obs.data <- obs.data[, c("latitude", "longitude")]
# Build species distribution model
bc.model <- bioclim(x = bioclim.data, p = obs.data)
## Error in bioclim(data.frame(m), ...): insufficient records
什么...？好的，此错误消息很难确定。但是，让我们考虑一下我们的obs.data数据框现在是什么样子：
head(obs.data)
## latitude longitude
## 1 32.33556 -110.8980
## 2 32.28267 -110.9028
## 3 30.29105 -110.7213
## 4 32.05413 -110.6837
## 5 32.25111 -110.7169
## 6 32.19404 -111.0198
第一列是纬度，第二列是经度，看起来不错。也就是说，直到我们考虑R通常如何处理坐标。当我们绘制东西时，通常使用如下语法：
plot(x, y)
需要注意的是，我们传递的第一个参数是x轴的数据，第二个参数是y轴的数据。该bioclim函数以相同顺序查找数据。也就是说，它将查看我们传递给的内容，p并假设第一列用于x轴，第二列用于y轴。但是我们的数据却是相反的顺序：第一列是纬度，本质上是y轴数据，第二列是经度，对应于x轴数据。因此，在传递obs.data给之前，我们需要反转列顺序bioclim：
# Reverse order of columns
obs.data <- obs.data[, c("longitude", "latitude")]
# Build species distribution model
bc.model <- bioclim(x = bioclim.data, p = obs.data)
呜呜！这里没有错误（希望如此）。
在地图上绘制模型之前，我们还需要采取进一步的步骤。我们需要生成一个对象，该对象的模型出现了柱仙人掌的概率。我们使用包装中的predict模型dismo：
# Predict presence from model
predict.presence <- dismo::predict(object = bc.model,
x = bioclim.data,
ext = geographic.extent)
您可能想知道为什么我们使用dismo::predict而不是仅仅使用predict。毫不奇怪，不同的程序包有时使用相同的函数名来执行非常不同的操作。在的情况下predict，有加载到内存中至少有三个包有一个predict函数：dismo，sp，和stats。尽管我们仅使用（R将使用版本）可能会很好，但是指定版本会将此事实显式地传达给其他阅读代码的人。因此，我们可以使用语法而不是让别人（或您将来的自己！）猜测。predictdismodismodismo::predict
足够！现在该作图了。我们像以前一样开始，使用空白的灰色地图添加模型，如果感觉喜欢，则将原始观测值添加为点。
# Plot base map
plot(wrld_simpl,
xlim = c(min.lon, max.lon),
ylim = c(min.lat, max.lat),
axes = TRUE,
col = "grey95")
# Add model probabilities
plot(predict.presence, add = TRUE)
# Redraw those country borders
plot(wrld_simpl, add = TRUE, border = "grey5")
# Add original observations
points(obs.data$longitude, obs.data$latitude, col = "olivedrab", pch = 20, cex = 0.75)
box()
该图显示了整个地图上出现柱仙人掌的可能性。请注意，所有值均低于1.0；实际上，根据该模型，地图上任何地方的最大概率仅为0.78。但是，我们很确定在索诺兰沙漠的相当宽广的地区都发现了柱仙人掌-毕竟，我们的观察结果证明了这一点！如果我们希望我们的地图更好地反映这一点，我们将需要重新运行分析，但是这次包括一些缺席点，据称这些地方不会出现柱仙人掌。问题是，我们只有柱仙人掌的存在数据。
伪缺席点
强制将仅在场数据与在场/不在场方法一起使用的一种常见解决方法是使用伪不在场或“背景”点。虽然“伪缺席”听起来很花哨，但这实际上仅意味着一个人从给定地理区域中随机采样点，并将其视为缺少感兴趣物种的位置。Barbet-Massin等人的研究是研究伪缺点的影响和最佳实践的重要资源。（2012）（有关完整详细信息，请参见下面的其他资源）。
为了我们的目的，我们将随机创建一组背景（又名伪缺）点，其数量与我们观察到的点数一样多。我们将使用bioclim数据文件确定点的空间分辨率，并将采样区域限制在柱仙人掌观测的一般区域。
# Use the bioclim data files for sampling resolution
bil.files <- list.files(path = "data/wc2-5",
pattern = "*.bil$",
full.names = TRUE)
# We only need one file, so use the first one in the list of .bil files
mask <- raster(bil.files[1])
# Randomly sample points (same number as our observed points)
background <- randomPoints(mask = mask, # Provides resolution of sampling points
n = nrow(obs.data), # Number of random points
ext = geographic.extent, # Spatially restricts sampling
extf = 1.25) # Expands sampling a little bit
快速浏览一下background我们刚刚创建的对象：
head(background)
## x y
## [1,] -108.6458 27.06250
## [2,] -109.6042 27.85417
## [3,] -112.2708 35.97917
## [4,] -109.1875 35.72917
## [5,] -111.9375 32.02083
## [6,] -112.8542 28.27083
我们还可以像在观察点上那样在地图上可视化它们：
# Plot the base map
plot(wrld_simpl,
xlim = c(min.lon, max.lon),
ylim = c(min.lat, max.lat),
axes = TRUE,
col = "grey95",
main = "Presence and pseudo-absence points")
# Add the background points
points(background, col = "grey30", pch = 1, cex = 0.75)
# Add the observations
points(x = obs.data$longitude,
y = obs.data$latitude,
col = "olivedrab",
pch = 20,
cex = 0.75)
box()
现在我们有了伪缺席点，我们需要再采取一步。要获得外观更传统的地图，需要对模型进行事后评估。为了进行评估，我们将仅使用部分数据（训练数据）构建模型，并保留一部分数据以在模型构建后评估（测试数据）。我们将保留20％的数据用于测试，因此我们使用包中的kfold函数dismo将每个观察值平均分配给一个随机组。
# Arbitrarily assign group 1 as the testing data group
testing.group <- 1
# Create vector of group memberships
group.presence <- kfold(x = obs.data, k = 5) # kfold is in dismo package
现在暂停一分钟，看看group.presence刚才创建的向量：
head(group.presence)
## [1] 1 1 4 2 3 5
# Should see even representation in each group
table(group.presence)
## group.presence
## 1 2 3 4 5
## 80 80 80 80 80
的输出table显示已将五个点分配给五个组中的每个组。在这种情况下，我们可以看到这些点已经均匀分布，其中第1组（测试组）中有20％的点。
我们将group.presence向量与观察到的数据结合使用，将观察到的数据分为训练数据集和测试数据集：
# Separate observations into training and testing groups
presence.train <- obs.data[group.presence != testing.group, ]
presence.test <- obs.data[group.presence == testing.group, ]
# Repeat the process for pseudo-absence points
group.background <- kfold(x = background, k = 5)
background.train <- background[group.background != testing.group, ]
background.test <- background[group.background == testing.group, ]
训练和测试模型
现在，我们有了（1）伪缺失点和（2）单独的训练和测试数据，我们可以重新构建模型，评估模型的性能并绘制出更加美观的地图。我们使用bioclim之前的功能构建模型，但是我们没有使用所有的观测值，obs.data而是仅使用存储在中的训练数据presence.train：
# Build a model using training data
bc.model <- bioclim(x = bioclim.data, p = presence.train)
# Predict presence from model (same as previously, but with the update model)
predict.presence <- dismo::predict(object = bc.model,
x = bioclim.data,
ext = geographic.extent)
现在，我们采用该模型，并使用观测数据和为模型测试保留的伪缺失点对其进行评估。然后，我们使用此测试确定发生概率的临界值，以确定柱仙人掌范围的边界。
# Use testing data for model evaluation
bc.eval <- evaluate(p = presence.test, # The presence testing data
a = background.test, # The absence testing data
model = bc.model, # The model we are evaluating
x = bioclim.data) # Climatic variables for use by model
# Determine minimum threshold for "presence"
bc.threshold <- threshold(x = bc.eval, stat = "spec_sens")
该threshold功能提供了多种通过stat参数确定阈值截止的方法。在这里，我们选择"spec_sens"，设置“灵敏度（真阳性率）和特异性（真阴性率）之和最高的阈值”。有关更多信息，请查看threshold（?threshold，还记得吗？）的文档。
而最后，我们可以把这个门槛绘制的地图与仙人掌的预测范围！
# Plot base map
plot(wrld_simpl,
xlim = c(min.lon, max.lon),
ylim = c(min.lat, max.lat),
axes = TRUE,
col = "grey95")
# Only plot areas where probability of occurrence is greater than the threshold
plot(predict.presence > bc.threshold,
add = TRUE,
legend = FALSE,
col = "olivedrab")
# And add those observations
points(x = obs.data$longitude,
y = obs.data$latitude,
col = "black",
pch = "+",
cex = 0.75)
# Redraw those country borders
plot(wrld_simpl, add = TRUE, border = "grey5")
box()
嗯...看起来不对。它绘制了地图绿色的很大一部分。让我们看一下我们实际要求R绘制的内容，即绘制的值predict.presence > bc.threshold。那是什么
predict.presence > bc.threshold
## class : RasterLayer
## dimensions : 216, 144, 31104 (nrow, ncol, ncell)
## resolution : 0.04166667, 0.04166667 (x, y)
## extent : -115, -109, 26, 35 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 0, 1 (min, max)
这两个栅格的比较将生成另一个值仅为0或1：0的栅格，其中比较的结果为FALSE（即，当栅格像元中的predict.presence值小于或等于对应栅格像元中的值时bc.threshold），并且1，比较结果为TRUE。由于此比较中有两个值（values字段中的0和1 ），因此我们需要更新col在图调用中传递给参数的值。我们不仅提供单个值，还提供了0（NA）和1（"olivedrab"）的颜色：
# Plot base map
plot(wrld_simpl,
xlim = c(min.lon, max.lon),
ylim = c(min.lat, max.lat),
axes = TRUE,
col = "grey95")
# Only plot areas where probability of occurrence is greater than the threshold
plot(predict.presence > bc.threshold,
add = TRUE,
legend = FALSE,
col = c(NA, "olivedrab"))
# And add those observations
points(x = obs.data$longitude,
y = obs.data$latitude,
col = "black",
pch = "+",
cex = 0.75)
# Redraw those country borders
plot(wrld_simpl, add = TRUE, border = "grey5")
box()
关于我们的方法的最后一点说明：我们绘制的地图对景观中的特定点是否适合感兴趣的物种进行了分类。这种分类在很大程度上取决于阈值（请参阅参考资料bc.threshold和文档threshold）和伪缺点。假设我们使用随机抽样生成了这些伪缺失点，那么如果您多次运行此代码，则预测范围内可能会出现变化（尝试！如果您从创建伪伪代码开始重新运行该代码，请尝试！缺席点，几乎可以保证您拥有一张不同的地图。）有很多方法可以解决这种变化，Barbet-Massin等人的论文也有此报道。（2012年）是很棒的资源 我将把它留作作业，供您确定哪种方法最合适！
我们的最终脚本，用于生成模型，确定阈值并可视化结果：
# Species distribution modeling for saguaro
# Jeff Oliver
# jcoliver@email.arizona.edu
# 2018-02-27
rm(list = ls())
# Load additional packages
library("sp")
library("raster")
library("maptools")
library("rgdal")
library("dismo")
# Download bioclim data and store in bioclim.data variable
bioclim.data <- getData(name = "worldclim",
var = "bio",
res = 2.5,
path = "data/")
# Read in saguaro observations
obs.data <- read.csv(file = "data/Carnegiea-gigantea-GBIF.csv")
# Drop any rows with NAs
obs.data <- obs.data[!is.na(obs.data$latitude), ]
# Only pull out those columns of interest and in the order we want them
obs.data <- obs.data[, c("longitude", "latitude")]
# Determine geographic extent of our data
max.lat = ceiling(max(obs.data$latitude))
min.lat = floor(min(obs.data$latitude))
max.lon = ceiling(max(obs.data$longitude))
min.lon = floor(min(obs.data$longitude))
geographic.extent <- extent(x = c(min.lon, max.lon, min.lat, max.lat))
# Crop the bioclim data to geographic extent of saguaro
bioclim.data <- crop(x = bioclim.data, y = geographic.extent)
# Create pseudo-absence, or background, points
# Use the bioclim data files for sampling resolution
bil.files <- list.files(path = "data/wc2-5",
pattern = "*.bil$",
full.names = TRUE)
# We only need one file, so use the first one in the list of .bil files
mask <- raster(bil.files[1])
# Randomly sample points (same number as our observed points)
background <- randomPoints(mask = mask, # Provides resolution of sampling points
n = nrow(obs.data), # Number of random points
ext = geographic.extent, # Spatially restricts sampling
extf = 1.25) # Expands sampling a little bit
# Arbitrarily assign group 1 as the testing data group
testing.group <- 1
# Create vector of group memberships
group.presence <- kfold(x = obs.data, k = 5) # kfold is in dismo package
# Separate observations into training and testing groups
presence.train <- obs.data[group.presence != testing.group, ]
presence.test <- obs.data[group.presence == testing.group, ]
# Repeat the process for pseudo-absence points
group.background <- kfold(x = background, k = 5)
background.train <- background[group.background != testing.group, ]
background.test <- background[group.background == testing.group, ]
# Build a model using training data
bc.model <- bioclim(x = bioclim.data, p = presence.train)
# Predict presence from model
predict.presence <- dismo::predict(object = bc.model,
x = bioclim.data,
ext = geographic.extent)
# Use testing data for model evaluation
bc.eval <- evaluate(p = presence.test, # The presence testing data
a = background.test, # The absence testing data
model = bc.model, # The model we are evaluating
x = bioclim.data) # Climatic variables for use by model
# Determine minimum threshold for "presence"
bc.threshold <- threshold(x = bc.eval, stat = "spec_sens")
# Load map data for plotting
data(wrld_simpl)
# Plot base map
plot(wrld_simpl,
xlim = c(min.lon, max.lon),
ylim = c(min.lat, max.lat),
axes = TRUE,
col = "grey95")
# Only plot areas where probability of occurrence is greater than the threshold
plot(predict.presence > bc.threshold,
add = TRUE,
legend = FALSE,
col = c(NA, "olivedrab"))
# And add those observations
points(x = obs.data$longitude,
y = obs.data$latitude,
col = "black",
pch = "+",
cex = 0.6)
# Redraw those country borders
plot(wrld_simpl, add = TRUE, border = "grey5")
box()
额外资源
dismo包装小插图
R中的物种分布模型
使用高斯过程的快速灵活的贝叶斯物种分布建模
运行一系列物种分布模型
Google地图上的SDM多边形
R包“ maxnet”，用于实现Java maxent包的功能
甲上伪缺勤在SDMS效果研究（Barbet的-Massin等人。2012）
一个PDF版本本课
返回学习者主页
有什么问题吗 给我发电子邮件jcoliver@email.arizona.edu。